#include "FEMSpace.h"
#include <typeinfo>
#include <math.h>
#include <ctime>

#define pi 4.0*atan(1.0)
double f(double *p)
{
    return 2 * pi*  pi* sin(pi * p[0]) * sin(pi * p[1]);
}

double bc(double *p)
{
    return sin(pi * p[0]) * sin(pi * p[1]);
}

double f2(double *p)
{
    return 2 * pi*  pi* cos(pi * p[0]) * cos(pi * p[1]);
    //return 1;
}

double bc2(double *p)
{
    return cos(pi * p[0]) * cos(pi * p[1]);
}
double a(double *p)
{
    return 1;
}
double g1(double *p)
{
    return exp(p[0] - 1.0);
    //return 0;
}
double g2(double *p)
{
    return exp(1.0 + p[1]);
    //return 0;
}
double g3(double *p)
{
    return exp(p[0] + 1.0);
    //return 0;
}
double g4(double *p)
{
    return exp(-1.0 + p[1]);
    //return 0;
}
double u(double *p)
{
    return exp(p[0] + p[1]);
    // double z = 0;
    // for(int k = 1; k < 200; k += 2)
    // {
	//     z += sin(k * pi * (1 + p[0]) / 2) / (k * k * k * sinh(k * pi)) * (sinh(k * pi * (1 + p[1]) / 2) + sinh(k * pi * (1 - p[1]) / 2));
    // }
    // return (1 - p[0] * p[0]) / 2 - 16 * z / (pi * pi * pi);
}
double fe(double *p)
{
    return -2 * exp(p[0] + p[1]);
    //return 1;
}
int main(int argc, char* argv[])
{
    //g++ -o main test_possion.cpp -std=c++11 -I /usr/include/eigen3/
    // RectangleDomain* r = new RectangleDomain({{-1,-1},{1,-1},{1,1},{-1,1}});
    // std::vector<Boundary<2> > B = r->boundary();
    vector<Point<2>> p = {{-1.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {-1.0, 1.0}};
    Rectangle2D domain;
    domain.setVertexList(p);
    vector<int> seg({std::atoi(argv[1]), std::atoi(argv[2])});
    Mesh<2>* m = new Q1Mesh(domain, seg);
    TemplateElement<2>* e = new Quadrilateral_1_Element();
    Equation<2>* equ = new PossionEquation<2>();
    equ -> SetRightHandsTermFunction(fe);
    equ -> SetBoundaryConditionFunction({g1, g2, g3, g4});
    // BoundaryFunction<2> * bf1 = new Dirichlet<2>(g1,B[0]);
    // BoundaryFunction<2> * bf2 = new Dirichlet<2>(g2,B[1]);
    // BoundaryFunction<2> * bf3 = new Dirichlet<2>(g3,B[2]);
    // BoundaryFunction<2> * bf4 = new Dirichlet<2>(g4,B[3]);
    // BoundaryCondition<2> bfc;
    // bfc.add(bf1);
    // bfc.add(bf2);
    // bfc.add(bf3);
    // bfc.add(bf4);
    Possion_2D possionproblem(m,e,equ);
    clock_t start,end;
    start = clock();
    possionproblem.setup_system();
    end=clock();		//程序结束用时
    double endtime=(double)(end-start)/CLOCKS_PER_SEC;
    cout<<"setup_system takes  "<<endtime*1000<<" ms"<<endl;	//ms为单位
    start = clock();
    possionproblem.assemble_system();
    end=clock();		//程序结束用时
    endtime=(double)(end-start)/CLOCKS_PER_SEC;
    cout<<"assemble_system takes  "<<endtime*1000<<" ms"<<endl;	//ms为单位
    start = clock();
    possionproblem.dealwithBoundaryCondition();
    end=clock();		//程序结束用时
    endtime=(double)(end-start)/CLOCKS_PER_SEC;
    cout<<"dealwithBoundaryCondition takes  "<<endtime*1000<<" ms"<<endl;	//ms为单位
    start = clock();
    possionproblem.solve();
    end=clock();		//程序结束用时
    endtime=(double)(end-start)/CLOCKS_PER_SEC;
    cout<<"solve takes  "<<endtime*1000<<" ms"<<endl;	//ms为单位
    // Eigen::VectorXd residual = possionproblem.A() * possionproblem.getSolution() - possionproblem.Rhs();
    // double residual_l2 = residual.squaredNorm();
    // double residual_max = residual.lpNorm<Eigen::Infinity>();
    // std::cout<<"l2 residual is :" << residual_l2 << " ,residual max norm is :" << residual_max<<std::endl;
    Eigen::VectorXd Truesolution = Eigen::MatrixXd::Zero(m->getTotalNumDofs(), 1);
    Point<2> p0 = domain.getVertex(0);
    Point<2> p2 = domain.getVertex(2);
    double xLen = (p2[0] - p0[0]) / (1.0 * seg[0]);
    double yLen = (p2[1] - p0[1]) / (1.0 * seg[1]);
    for(int j = 0; j < seg[1] + 1; ++j)
    {
	for(int i = 0; i < seg[0] + 1; ++i)
	{
	    //s[i+j*(segmentx+1)] = (TrueSolve(p0[0]+i*xLen,p0[0]+j*yLen));//stupid!!!
	    double point[2] = {p0[0] + i * xLen, p0[0] + j * yLen};
	    Truesolution[i + j * (seg[0] + 1)] = u(point);//stupid!!!
	}
    }
    // for(int j = 0; j <= seg[1]; j++)
    // {
    //     for(int i = 0;i <= seg[0]; i++)
    //     {
    //         double point[2] = {-1 + i * m->x_h(),-1 + j * m->y_h()};
    //         Truesolution[i + (seg[0] + 1) * j] = u(point);
    //     }
    // }
    // std::cout << "done" << std::endl;
    Eigen::VectorXd error = Truesolution - possionproblem.getSolution();
    //std::cout << "error is :" << error <<std::endl;
    double error_l2 = error.norm();
    double error_max = error.lpNorm<Eigen::Infinity>();
    std::cout<<"Error l2 norm is :" << error_l2 << ", error max norm is :" << error_max<<std::endl;
    //
    if(std::atoi(argv[3]) == 0)
    {
        FILE* file1;
        file1 = fopen("output.txt", "w");//“w”写文件覆盖原内容，不想覆盖则用“a”
        if (NULL == file1)
    	    printf("error");
        fprintf(file1, "%f\n%f\n%f\n%f\n", domain.getVertex(0)[0], domain.getVertex(1)[0], domain.getVertex(0)[1], domain.getVertex(3)[1]);
        fprintf(file1, "%d\n%d\n%d\n", m->getSegment()[0], m->getSegment()[1], 1);
        for (int i = 0; i < possionproblem.getSolution().size(); i++)
    	    fprintf(file1,"%.12f\n", possionproblem.getSolution()[i]);//循环输出数据并写入
        fclose(file1);//最后一定记得关闭文件
        std::cout << "Data has been output" << std::endl; 
    }
    else if(std::atoi(argv[3]) == 1)
    {
        FILE* file2;
        file2 = fopen("err.csv", "a");
        if (NULL == file2)
    	    printf("error");
        fprintf(file2, "%d %d %.12f\n", std::atoi(argv[1]), std::atoi(argv[2]), error_l2);
        std::cout << "Data has been output" << std::endl; 
    }
    else
        std::cout << "No output" << endl;
}
